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1. Introduction 

The anomalous magnetic moment of the muon is one of the most precisely measured quantities 
in particle physics and it provides a stringent test of the standard model. The current experimental 
measurement and theoretical prediction of the muon anomalous magnetic moment = (g —2)^/2 
show a tension ranging between 3 and 4 standard deviations [1,2]. With the planned improvement 
of the experimental precision at Fermilab and at J-PARC, further reduction of the theoretical uncer¬ 
tainty is required, in order to be able to resolve the current discrepancy and potentially gain insight 
into new physics. The dominant uncertainties in the theoretical prediction of are of hadronic ori¬ 
gin: the error of the leading hadronic contribution and the error of the hadronic light by light 

contribution Since it was worked out that the hadronic vacuum polarization (HVP) can be 

evaluated in Euclidean space-time [3], the extraction of from first principles has become one 
of the long-standing goals of lattice community (for a recent review see [4]). 

The leading hadronic contribution to is obtained as 

C° = (^)'r^2W)xn(e2), (1.1) 

^ n Jo 

where = n(2^) — n(0) denotes the renormalized vacuum polarization tensor and f{Q^) is 

a known function diverging as —)• 0, so that the integrand is strongly peaked at « m^/4. 
Although the HVP function is accessible from lattice simulations, the total integral is dominated in 
the low momentum region where relative errors of Tl{Q^) are enhanced. On top of that, the HVP 
determinations from the lattice show significant pion mass dependence [10], as well as significant 
decrease in the signal to noise ratio near 2^ = 0 as —)• All this makes the control the 

systematics in lattice determination of challenging even before performing the continuum 
extrapolation and including the finite volume effects, isospin breaking effects etc. 

A traditional approach to extracting from lattice computations involves assuming a func¬ 
tional dependence for Yl{Q^) and performing a fit of the lattice data to obtain the infrared subtrac¬ 
tion n(0). Let us note that a number of recently proposed methods allow for direct computation 
of n(0) [11, 12], where the method of Ref. [11] requires additional computational effort that was 
not affordable with DWF fermions at the physical point and the application of [12] to the DWF 
ensembles still requires further investigation. In this work, we stick to the standard approach which 
requires the extrapolation to n(0). We eliminate the systematics from chiral extrapolation by mea¬ 
suring HVP directly at the physical quark mass and apply a series of Pade approximants [13] in an 
attempt to asses the systematics introduced by fitting. 

2. DWF action, algorithms and computational strategy 

In this preliminary work we use a single ensemble of 2 -|- 1 flavors of Mobius domain wall 
fermions [5] wifh physical pion mass and lattice extent 48^ x 96 x 24 (see the left panel of Figure 
1). The ensemble generation has been done with the CPS QCD* package as a part of the wide 
spread effort of the RBC-UKQCD collaboration. We quote here only the total number of measured 
configurations Nconf = 88 separated by 20 molecular dynamics trajectories, the lattice spacing for 

'Columbia Physics Systems, http://qcdoc.phys.columbia.edu/cps.html 
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this ensemble = 1.7295(38)GeV and the renormalization constant for the vector current Zy = 
0.71076(25); for more details on the RBC-UKQCD DWF ensembles see [8]. The measurements 
of the HVP tensor are using the BAGEL^ assembler kernel generator and the UKHadron codes. 
Its recent adaptation includes the hierarchically deflated conjugate gradient (HDCG) algorithm [6]. 
The usage of the HDCG algorithm on this particular ensemble gave us a 3.5 x speed up over EigCG 
and made this computation possible with the available computer resources due to the vastly reduced 
memory footprint (64 vs. 1500 vectors). 

Our strategy for the extraction of the HVP form factor involved the computation of the local- 
conserved vector two point function, as previously done in [9, 10] 

( 2 . 1 ) 

where denotes conserved Mobius current, Jy denotes the local current and Q are lattice mo¬ 
menta^. The conserved current for the Mobius DWE action and the details of its implementation 
are given in [7, 8]. The choice of local-conserved (E-C) correlator leaves one Ward identity intact 
and it is favourable over the computation of the local-local (E-E) and conserved-conserved (C-C) 
correlator. Namely, C-E still leaves the Ward identity QpTlfiv{Q) = 0 intact, while E-E on the 
lattice does not satisfy a WI for any of the indices, C-C has an additional contact term that needs 
to be subtracted and in addition C-C is computationally more expensive than the E-C construction 
for DWE. 

We can then relate the HVP form factor with the HVP tensor 

n^v{Q) = {5pvQ^-QpQvMQ^). ( 2 . 2 ) 

We work with the diagonal components of the HVP tensor, n^^(2), and furthermore only inlcude 
in the analysis momenta which have no component in the /r direction, = 0. This selection 
simplifies the calculation of the HVP form factor ^ eliminates a fraction of the 

unwanted cutoff effects [9]. 

3. Point source vs. stochastic sources 

Extended stochastic sources have often been proven to be advantageous in the computation 
of the two point correlators [14], giving the significantly better statistical precision than the point 
source for same computational cost. The ’one-end trick’ introduced in [14, 15] and tested for DWE 
in [16] has been a standard choice of RBC-UKQCD for the calculation of meson two-point and 
three-point functions. In order to find an optimal setup for the HVP computation on our physical 
point ensemble of interest, we carry out the same cost comparison of its computation with a Z(2) x 
Z(2) stochastic wall source [16] with a more conventional computation of the HVP with a point 
source. We find that after performing the Eourier transform in (2.1), the volume source performs 
significantly better than the point source in the large momentum region, but below momenta k, 
0.2GeV^ the precision of the point outperforms the volume source we tested. This is illustrated in 

^http://www. ph.ed.ac.uk/ paboyle/bagel/Bagel.html 

^The components of the lattice momenta Q are defined as Q^i = 
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Figure 1: Left: overview of the = 2 + 1 RBC-UKQCD ensembles used in the HVP programme. Full red 
circle denotes an ensemble used in this work, its lattice size being 48^ x 96 x 24 and m^ = 138GeV. Right: 
comparison of the relative precision of the scalar vacuum polarisation computed with a point source (squares) 
and a Z(2) X Z(2) stochastic source (triangles) averaged over 12 timeslices on a 48^ x 96 x 24 physical point 
DWF ensemble. On the y-axis we plot values of 511(2) = An(2)/n(2) obtained by the same number of 
source positions with both source types, while the x-axis reads lattice momenta in GeV^. 


the right panel of Figure 1, where the relative statistical precision of Tl{Q^) for the two choices of 
the source averaged over 12 timeslices is shown. Since this low momentum region is the one where 
rougly 90% of the contribution to comes from [19], we expect the precision of the lattice 
data in this momentum range to be the most important for achieving the desired precision on 
We therefore use point source data for the analysis presented in the following. It is possible that 
the noise characteristics may differ if different analysis methods are used [17] and that will be the 
subject of further study. 


4. Fade fits 


The evaluation of the integral (1.1) requires asuming a functional dependence of 11(2^). In 
the first attempts to obtain from the lattice, many of the functional forms used to date were 
based on Vector Meson Dominance (VMD). Recent studies, however, provide strong evidence that 
this approach introduces uncontrolled systematics [13, 18]. We use in this work a series of Fade 
fits 


n[A?,D](2^) 




(4.1) 


to test the stability of the fit for different number of fit parameters, without relying on VMD. We 
use the physical pion mass data set described in Section 2 and perform [1,1], [2,1] and [2,2] Fade 
fits which have 3, 4 and 5 free parameters respectively. These lie in the sequence of Fades proposed 
in Ref. [13], known to be such that the full sequence converges to the actual polarization for any 
compact region in the complex plane excluding the cut along the negative real axis. 

Let us first discuss the extrapolation of to = 0. In Figure 2 we show different orders 

of Fades applied to a fitting interval [0, Q^]. For the lowest cut value we took = 1.5GeV^, the 
fit [2,2] gives large uncertainities and the extrapolated values of 11(0) from all attempted fits are 
compatible with each other. For the remaining three values of both Fades [2,1] and [2,2] give 
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Figure 2: Using different Fades for extrapolating 11(2^) to = 0 on a 48^ x 96 x 24 lattice with = 
1.7295GeV and m^ — 0.138GeV. Upper left plot includes in the fit momenta up to 2^ < 2c ~ l -5GeV^, 
upper right includes 2^ < 2c = 2GeV^, lower left plot includes 2^ < 2c = 2.5GeV^ in the fits and lower 
right fits momenta up to 2^ < 2c = 3GeV^. 


good x^-: while Fade [1,1] gives n(0) more than one sigma away from the higher order Fades. We 
ilustrate this in the remaining three panels of Figure 2. For cut values 2c = 2GeV^;2.5GeV^ and 
3GeV^ the extrapolated n(0) values from [2,1] and [2,2] Fades are compatible with each other and 
one might be tempted to consider [2,1] a fit of choice, since it has less free parameters then [2,2] 
and therefore gives better statistical precision. We will see in the following section that using lower 
order Fades in combination with higher cuts in the momenta 2c introduces systematics one would 
like to avoid. 

5. Including large in the Fade fits 

Including momenta up to higher values of 2c Pade fits with small number of parameters is 
yet another potential source of systematic uncertainty. This was first observed for the dispersive 
model constructed for studying the uncertainities associated with different fits to lattice data for the 
vacuum polarization [18]. Namely, when the maximal momenta included in the fits of the 7=1 
polarization data has been changed from 2c = IGeV^ to 2c = 1.5GeV^, the value of the relevant 
vacuum polarization contribution is undershooting its true value, which in this particular model 
can be directly extracted from the hadronic T-decays data (see right panel of Figure 3). We probe 
on our data set several values of 2c together with the Fade fits used in the previous section. In 
the left panel of Figure 3 we show the contributions < IGeV^) resulting from different 

Fade fits on our data set. Only higher order Fades seem to give consistent contributions to 
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Figure 3: Left: contribution to when the corresponding Fade fits defined in eq. (4.1) are performed in 
the momenta region G [0,2^] for different values of maximal included in the fit. Right: Result for 
the 7=1 non-Strange polarization function obtained from the dispersive model defined in [18]. The plotted 
contributions to 5^ from the momenta region [0, l]GeV^ for different Fade approximants^ taken 

from Table 1 in Ref. [18]. 


for higher values of as already suggested by the dispersive model study. In particular, on our 
ensemble and with the current statistics, both [1,1] and [2,1] Fade are not acceptable for the values 
of Qc ^ l-5GeV^. Even though the values of < IGeV^) obtained by [2,1] Fade function 

look compatible for different values of note that the plotted values are correlated and once the 
correlations are subtracted, the spread comes out to be larger than the corresponding estimates of 
the statistical errors which favours either using higher order Fades such as [2,2] or some kind of 
hybrid methods along the lines of [19] where the fitting interval is significantly shorter. 


6. Summary 

We are measuring the leading order hadronic contribution to with 2 + 1 flavors of DWF 
fermions af fhe physical poinf. Af fhis poinf of fhe calculation and on fhe sfudied physical poinf dafa 
sef, we find fhaf fhe poinf source gives beffer relafive precision fhan fhe sfochasfic source wifh fhe 
one-end frick for fhe momenfa < 0.2GeV^, which is fhe region fhaf gives mosf of fhe fofal 
confribufion. We perform fifs using Fade approximanfs, which consfifufe a model-independenf 
paramefrisafion of fhe HVF form factor, and fesf how changing fhe number of free paramefers 
and fhe tiffing range affecfs fhe infrared subfracfion and fhe evaluafion of fhe infegral (1.1). The 
Fades we employ lie in fhe sequence proposed in Ref. [13] and are guaranfeed to converge fo 
fhe acfual polarization, buf fhe sysfemafics associafed wifh fhe use of parficular order Fades wifh 
particular choices of remain to be explored in more defail. Frogress sfabilizing fhe fif resulfs 
may also be possible using fhe conformal variable polynomial fif forms discussed in Ref. [18]. 
The confribufions fo we are gelling in Seclion 5 are in fhe same ballpark as previous HVF 
compulalions, buf wifh large sfatisfical uncerfainilies. Therefore, already wifh fhis compufafion on 
a single ensemble we can see fhaf more sophislicaled mefhods are necessary fo exlracl fhe value 
wifh fhe precision comparable fo fhe experimenlal one. 

^Note that Fade approximants [A, 77] plotted in the right panel of Figure 3 correspond to [N-1,D] in the notation of 
Ref. [18]. 
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